Load packages

### 0. Load packages 
if(!require(qgraph)) install.packages(pkgs="qgraph",repos="http://cran.r-project.org")

## Loading required package: qgraph

## Registered S3 methods overwritten by 'huge':
##   method    from   
##   plot.sim  BDgraph
##   print.sim BDgraph

if(!require(bootnet)) install.packages(pkgs="bootnet",repos="http://cran.r-project.org")

## Loading required package: bootnet

## Loading required package: ggplot2

## This is bootnet 1.4.3

## For questions and issues, please see github.com/SachaEpskamp/bootnet.

if(!require(dplyr)) install.packages(pkgs="dplyr",repos="http://cran.r-project.org")

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

if(!require(igraph)) install.packages(pkgs="igraph",repos="http://cran.r-project.org")

## Loading required package: igraph

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

if(!require(doParallel)) install.packages(pkgs="doParallel",repos="http://cran.r-project.org")

## Loading required package: doParallel

## Loading required package: foreach

## Loading required package: iterators

## Loading required package: parallel

library("foreign")
library("qgraph")
library("bootnet")
library("dplyr")
library("igraph")
library("doParallel")

Generate true network

### 1. Model selection example (figure 1 & 2)

# Set random seed:
set.seed(3)

# Generate true network:
# relationship between nodes
# an adjacency matrix: a square matrix used to represent a finite graph
Kappa <- as.matrix(get.adjacency(watts.strogatz.game(1,8,1,0)))
Kappa 

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    0    1    0    0    0    0    0    1
## [2,]    1    0    1    0    0    0    0    0
## [3,]    0    1    0    1    0    0    0    0
## [4,]    0    0    1    0    1    0    0    0
## [5,]    0    0    0    1    0    1    0    0
## [6,]    0    0    0    0    1    0    1    0
## [7,]    0    0    0    0    0    1    0    1
## [8,]    1    0    0    0    0    0    1    0

help("watts.strogatz.game")

The Watts-Strogatz small-world model
Description
Generate a graph according to the Watts-Strogatz network model.

Usage
sample_smallworld(dim, size, nei, p, loops = FALSE, multiple = FALSE)

smallworld(...)
Arguments
dim 
Integer constant, the dimension of the starting lattice.

size    
Integer constant, the size of the lattice along each dimension.

nei 
Integer constant, the neighborhood within which the vertices of the lattice will be connected.

p   
Real constant between zero and one, the rewiring probability.

as.matrix(get.adjacency(sample_smallworld(1,8,1,0)))

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    0    1    0    0    0    0    0    1
## [2,]    1    0    1    0    0    0    0    0
## [3,]    0    1    0    1    0    0    0    0
## [4,]    0    0    1    0    1    0    0    0
## [5,]    0    0    0    1    0    1    0    0
## [6,]    0    0    0    0    1    0    1    0
## [7,]    0    0    0    0    0    1    0    1
## [8,]    1    0    0    0    0    0    1    0

# coefficients between nodes (upper.tri)
Kappa[upper.tri(Kappa)] <- runif(sum(upper.tri(Kappa)),0.3,0.4) * Kappa[upper.tri(Kappa)]  
# coefficients between nodes (upper.tri & lower.tri)
Kappa[lower.tri(Kappa)] <- t(Kappa)[lower.tri(Kappa)] 
Kappa

##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## [1,] 0.0000000 0.3168042 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] 0.3168042 0.0000000 0.3384942 0.0000000 0.0000000 0.0000000 0.0000000
## [3,] 0.0000000 0.3384942 0.0000000 0.3604394 0.0000000 0.0000000 0.0000000
## [4,] 0.0000000 0.0000000 0.3604394 0.0000000 0.3630979 0.0000000 0.0000000
## [5,] 0.0000000 0.0000000 0.0000000 0.3630979 0.0000000 0.3867919 0.0000000
## [6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.3867919 0.0000000 0.3228202
## [7,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3228202 0.0000000
## [8,] 0.3015330 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3910148
##           [,8]
## [1,] 0.3015330
## [2,] 0.0000000
## [3,] 0.0000000
## [4,] 0.0000000
## [5,] 0.0000000
## [6,] 0.0000000
## [7,] 0.3910148
## [8,] 0.0000000

# diagnoal parameters
Kappa <- -1*Kappa
diag(Kappa) <- 1
Kappa

##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,]  1.0000000 -0.3168042  0.0000000  0.0000000  0.0000000  0.0000000
## [2,] -0.3168042  1.0000000 -0.3384942  0.0000000  0.0000000  0.0000000
## [3,]  0.0000000 -0.3384942  1.0000000 -0.3604394  0.0000000  0.0000000
## [4,]  0.0000000  0.0000000 -0.3604394  1.0000000 -0.3630979  0.0000000
## [5,]  0.0000000  0.0000000  0.0000000 -0.3630979  1.0000000 -0.3867919
## [6,]  0.0000000  0.0000000  0.0000000  0.0000000 -0.3867919  1.0000000
## [7,]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 -0.3228202
## [8,] -0.3015330  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
##            [,7]       [,8]
## [1,]  0.0000000 -0.3015330
## [2,]  0.0000000  0.0000000
## [3,]  0.0000000  0.0000000
## [4,]  0.0000000  0.0000000
## [5,]  0.0000000  0.0000000
## [6,] -0.3228202  0.0000000
## [7,]  1.0000000 -0.3910148
## [8,] -0.3910148  1.0000000

# Make 4 and 5 negative
Kappa[1:3,4:5] <- -1 * Kappa[1:3,4:5]
Kappa[4:5,1:3] <- -1 * Kappa[4:5,1:3]

# True cov structure:
Sigma <- solve(Kappa)


# Sample data:
library("mvtnorm")
Data <- rmvnorm(100,sigma = Sigma)

Model estimation

# Correlation matrix:
corMat <- cor(Data)
corMat

##             [,1]       [,2]        [,3]        [,4]        [,5]        [,6]
## [1,]  1.00000000  0.5123636  0.14843729 -0.06473797 -0.03746480 -0.03790132
## [2,]  0.51236364  1.0000000  0.37934557 -0.23037365 -0.11652514  0.11409230
## [3,]  0.14843729  0.3793456  1.00000000 -0.42766411 -0.19930338 -0.06559148
## [4,] -0.06473797 -0.2303736 -0.42766411  1.00000000  0.52231864  0.21352558
## [5,] -0.03746480 -0.1165251 -0.19930338  0.52231864  1.00000000  0.33790299
## [6,] -0.03790132  0.1140923 -0.06559148  0.21352558  0.33790299  1.00000000
## [7,]  0.17674246  0.1262173 -0.06112605  0.19761716  0.20360528  0.36757685
## [8,]  0.48707457  0.2431745  0.03451387  0.15583973  0.06445881  0.06435186
##             [,7]       [,8]
## [1,]  0.17674246 0.48707457
## [2,]  0.12621730 0.24317454
## [3,] -0.06112605 0.03451387
## [4,]  0.19761716 0.15583973
## [5,]  0.20360528 0.06445881
## [6,]  0.36757685 0.06435186
## [7,]  1.00000000 0.39751551
## [8,]  0.39751551 1.00000000

# Labels:
Labs <- LETTERS[1:ncol(Kappa)]

### Figure 1: True network structure
# pdf("Fig1.pdf",width=4,height=4)
qgraph(qgraph:::wi2net(Kappa), labels = Labs, vsize = 15, edge.labels=TRUE, vTrans = 254, edge.label.cex=1.5, theme = "colorblind")

# dev.off()


# EBIC analysis:
res0 <- EBICglasso(corMat,nrow(Data),0,returnAllResults=TRUE,nlambda=10)
res0.25 <- EBICglasso(corMat,nrow(Data),0.25,returnAllResults=TRUE,nlambda=10)
res0.5 <- EBICglasso(corMat,nrow(Data),0.5,returnAllResults=TRUE,nlambda=10)

### Figure 2: Estimated networks:
# pdf("Fig2.pdf",width=5*2,height=2*2)
layout(rbind(1:5,6:10))
for (i in 1:10){
  qgraph(qgraph::wi2net(res0$results$wi[,,i]), labels = Labs, vsize = 15,mar=c(2,2,2,2), theme = "colorblind")
  text(0,0.3,paste0("Lambda: ",round(res0$lambda[i],3)),cex=1)
  text(0,0.1,paste0("EBIC (gamma = 0): ",round(res0$ebic[i],1)),cex=ifelse(res0$ebic[i] == min(res0$ebic),1,0.6),font = ifelse(res0$ebic[i] == min(res0$ebic),2,1))
  text(0,-0.1,paste0("EBIC (gamma = 0.25): ",round(res0.25$ebic[i],1)),cex=ifelse(res0.25$ebic[i] == min(res0.25$ebic),1,0.6),font = ifelse(res0.25$ebic[i] == min(res0.25$ebic),2,1))
  text(0,-0.3,paste0("EBIC (gamma = 0.5): ",round(res0.5$ebic[i],1)),cex=ifelse(res0.5$ebic[i] == min(res0.5$ebic),1,0.6),font = ifelse(res0.5$ebic[i] == min(res0.5$ebic),2,1))
}

# dev.off()

Conditional Element Selection
Description
ifelse returns a value with the same shape as test which is filled with elements selected from either yes or no depending on whether the element of test is TRUE or FALSE.

Usage
ifelse(test, yes, no)
Arguments
test    
an object which can be coerced to logical mode.

yes 
return values for true elements of test.

no  
return values for false elements of test.

Model estimation & plot networks

##### 3 ESTIMATE & PLOT NETWORKS (Figure 3)

graph1 <- estimateNetwork(
  data,
  default = "EBICglasso",
  corMethod = "cor_auto",
  tuning = 0)

## Estimating Network. Using package::function:
##   - qgraph::EBICglasso for EBIC model selection
##     - using glasso::glasso
##   - qgraph::cor_auto for correlation computation
##     - using lavaan::lavCor

## Variables detected as ordinal: 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20

## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret the
## presence of the smallest edges with care. Setting threshold = TRUE will enforce
## higher specificity, at the cost of sensitivity.

#default ="EBICglasso" specifies that the EBICglasso function from qgraph is used

graph2 <- estimateNetwork(
  data,
  default = "EBICglasso",
  corMethod = "cor_auto",
  tuning = 0.25)

## Estimating Network. Using package::function:
##   - qgraph::EBICglasso for EBIC model selection
##     - using glasso::glasso
##   - qgraph::cor_auto for correlation computation
##     - using lavaan::lavCor
## Variables detected as ordinal: 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20

## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret the
## presence of the smallest edges with care. Setting threshold = TRUE will enforce
## higher specificity, at the cost of sensitivity.

graph3 <- estimateNetwork(
  data,
  default = "EBICglasso",
  corMethod = "cor_auto",
  tuning = 0.5)

## Estimating Network. Using package::function:
##   - qgraph::EBICglasso for EBIC model selection
##     - using glasso::glasso
##   - qgraph::cor_auto for correlation computation
##     - using lavaan::lavCor
## Variables detected as ordinal: 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20

#The corMethod = "cor_auto" argument specifies polychoric and polyserial correlations

# Compute average layout:
# To compare two networks, one should constrain the layout to be equal for both networks. One way to do so is by using averageLayout
L <- averageLayout(graph1, graph2, graph3)

#pdf("Fig3.pdf")
# maximum = 0.54 so it is consistent with later graphs 
plot(graph1, layout=L, maximum = 0.54, cut = 0) 

plot(graph2, layout=L, maximum = 0.54, cut = 0)

plot(graph3, layout=L, maximum = 0.54, cut = 0)

#dev.off()

wm1 <- getWmat(graph1)
wm2 <- getWmat(graph2)
wm3 <- getWmat(graph3)
wm1

##               1           2          3           4          5           6
## 1   0.000000000  0.18874500 0.02041634  0.26949321 0.09830280  0.00000000
## 2   0.188745000  0.00000000 0.41772177  0.00000000 0.02984737  0.00000000
## 3   0.020416337  0.41772177 0.00000000  0.21831074 0.13159891  0.00000000
## 4   0.269493205  0.00000000 0.21831074  0.00000000 0.28643618  0.04458847
## 5   0.098302796  0.02984737 0.13159891  0.28643618 0.00000000  0.18003268
## 6   0.000000000  0.00000000 0.00000000  0.04458847 0.18003268  0.00000000
## 7   0.000000000  0.06858381 0.02724846  0.12789258 0.12053477  0.21212668
## 8   0.000000000  0.00000000 0.00000000  0.00000000 0.00000000  0.09115236
## 9   0.000000000 -0.12634394 0.00000000  0.03406929 0.00000000 -0.08216170
## 10  0.000000000  0.00000000 0.00000000  0.09573875 0.00654187  0.03321589
## 11  0.194457214  0.03799185 0.02511442  0.06800126 0.01263520  0.00000000
## 12  0.062786336  0.03125524 0.09709947 -0.03650199 0.00000000  0.00000000
## 13  0.000000000  0.00000000 0.00000000 -0.03771817 0.00000000  0.00000000
## 14  0.000000000  0.00000000 0.00000000 -0.10174362 0.00000000  0.00000000
## 15  0.002444043  0.00000000 0.00000000  0.00000000 0.11178363  0.00000000
## 16 -0.004295477  0.14300460 0.00000000 -0.04662989 0.00000000  0.00000000
## 17  0.000000000  0.00000000 0.07662094  0.00000000 0.00000000  0.06282331
## 18 -0.026095863  0.07743453 0.08133495  0.00000000 0.11272345  0.00000000
## 19  0.025083092  0.00000000 0.00000000 -0.01111239 0.00000000  0.00000000
## 20  0.005666384  0.11676674 0.09682865  0.00000000 0.03893225  0.00000000
##              7          8           9          10         11          12
## 1   0.00000000 0.00000000  0.00000000  0.00000000 0.19445721  0.06278634
## 2   0.06858381 0.00000000 -0.12634394  0.00000000 0.03799185  0.03125524
## 3   0.02724846 0.00000000  0.00000000  0.00000000 0.02511442  0.09709947
## 4   0.12789258 0.00000000  0.03406929  0.09573875 0.06800126 -0.03650199
## 5   0.12053477 0.00000000  0.00000000  0.00654187 0.01263520  0.00000000
## 6   0.21212668 0.09115236 -0.08216170  0.03321589 0.00000000  0.00000000
## 7   0.00000000 0.03972146  0.00000000  0.00000000 0.00000000  0.14892145
## 8   0.03972146 0.00000000  0.01061640  0.08314477 0.02193945  0.00000000
## 9   0.00000000 0.01061640  0.00000000  0.01903501 0.27003465  0.17438501
## 10  0.00000000 0.08314477  0.01903501  0.00000000 0.47855705 -0.03917531
## 11  0.00000000 0.02193945  0.27003465  0.47855705 0.00000000  0.00000000
## 12  0.14892145 0.00000000  0.17438501 -0.03917531 0.00000000  0.00000000
## 13  0.00000000 0.00000000  0.04846113  0.00000000 0.01398274  0.13173568
## 14  0.03207787 0.00000000  0.13011228 -0.02027354 0.00000000  0.09084672
## 15  0.08373357 0.01699078  0.10368434  0.00000000 0.16221888  0.06550974
## 16  0.06201111 0.11061140  0.00000000  0.05790148 0.00000000  0.00000000
## 17  0.00000000 0.05797455  0.05298531  0.00000000 0.00000000  0.00000000
## 18  0.00000000 0.00000000  0.00000000 -0.03159653 0.00000000  0.01513076
## 19 -0.11817257 0.10797522  0.05759475  0.01598574 0.02382759  0.30016698
## 20  0.00000000 0.00000000  0.00000000  0.00000000 0.00000000  0.00000000
##             13          14          15           16         17          18
## 1   0.00000000  0.00000000 0.002444043 -0.004295477 0.00000000 -0.02609586
## 2   0.00000000  0.00000000 0.000000000  0.143004603 0.00000000  0.07743453
## 3   0.00000000  0.00000000 0.000000000  0.000000000 0.07662094  0.08133495
## 4  -0.03771817 -0.10174362 0.000000000 -0.046629892 0.00000000  0.00000000
## 5   0.00000000  0.00000000 0.111783627  0.000000000 0.00000000  0.11272345
## 6   0.00000000  0.00000000 0.000000000  0.000000000 0.06282331  0.00000000
## 7   0.00000000  0.03207787 0.083733569  0.062011109 0.00000000  0.00000000
## 8   0.00000000  0.00000000 0.016990779  0.110611398 0.05797455  0.00000000
## 9   0.04846113  0.13011228 0.103684338  0.000000000 0.05298531  0.00000000
## 10  0.00000000 -0.02027354 0.000000000  0.057901477 0.00000000 -0.03159653
## 11  0.01398274  0.00000000 0.162218881  0.000000000 0.00000000  0.00000000
## 12  0.13173568  0.09084672 0.065509741  0.000000000 0.00000000  0.01513076
## 13  0.00000000  0.47295317 0.025273362  0.061837398 0.08214261  0.09083379
## 14  0.47295317  0.00000000 0.113112038  0.000000000 0.00000000  0.00000000
## 15  0.02527336  0.11311204 0.000000000  0.266780501 0.00000000  0.00000000
## 16  0.06183740  0.00000000 0.266780501  0.000000000 0.03665648  0.11716494
## 17  0.08214261  0.00000000 0.000000000  0.036656477 0.00000000  0.32326412
## 18  0.09083379  0.00000000 0.000000000  0.117164939 0.32326412  0.00000000
## 19  0.15079596  0.02959822 0.012002099  0.000000000 0.00000000  0.09671241
## 20  0.11331185  0.02338779 0.000000000  0.000000000 0.01003874  0.10200155
##             19          20
## 1   0.02508309 0.005666384
## 2   0.00000000 0.116766744
## 3   0.00000000 0.096828655
## 4  -0.01111239 0.000000000
## 5   0.00000000 0.038932251
## 6   0.00000000 0.000000000
## 7  -0.11817257 0.000000000
## 8   0.10797522 0.000000000
## 9   0.05759475 0.000000000
## 10  0.01598574 0.000000000
## 11  0.02382759 0.000000000
## 12  0.30016698 0.000000000
## 13  0.15079596 0.113311854
## 14  0.02959822 0.023387794
## 15  0.01200210 0.000000000
## 16  0.00000000 0.000000000
## 17  0.00000000 0.010038744
## 18  0.09671241 0.102001553
## 19  0.00000000 0.201302918
## 20  0.20130292 0.000000000

sum(wm1[upper.tri(wm1, diag=F)]!=0)  #105 edges (85 missing)

## [1] 105

sum(wm2[upper.tri(wm2, diag=F)]!=0)  #95 edges (95 missing)

## [1] 95

sum(wm3[upper.tri(wm3, diag=F)]!=0)  #87 edges (103 missing)

## [1] 87

sum(abs(wm1[upper.tri(wm1, diag=F)])) /2 #5.091125 absolute sum 

## [1] 5.091125

sum(abs(wm2[upper.tri(wm2, diag=F)])) /2 #4.625315

## [1] 4.625315

sum(abs(wm3[upper.tri(wm3, diag=F)])) /2 #4.337963

## [1] 4.337963

##### 4 CENTRALITY (Figure 4)

#pdf("Fig4.pdf")
centralityPlot(list(EBIC0=graph1,EBIC0.25=graph2,EBIC0.5=graph3), include="all")

## Note: z-scores are shown on x-axis rather than raw centrality indices.

#dev.off()

# node strength, quantifying how well a node is directly connected to other nodes (该节点直接相连的连线加权值的绝对值之和)
# closeness, quantifying how well a node is indirectly connected to other nodes (所有节点到该节点的最短距离之和的倒数)
# betweenness, quantifying how important a node is in the average path between two other nodes. (该节点出现在网络中任意两个节点的最短路径上的次数)

Prior sample size analysis

trueNetwork <- estimateNetwork(
  data,
  default = "EBICglasso",
  corMethod = "cor_auto",
  tuning = 0.5,
  refit = TRUE)

## Estimating Network. Using package::function:
##   - qgraph::EBICglasso for EBIC model selection
##     - using glasso::glasso
##   - qgraph::cor_auto for correlation computation
##     - using lavaan::lavCor

## Variables detected as ordinal: 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20

## Refitting network without LASSO regularization

library("bootnet")
simRes <- netSimulator(trueNetwork$graph, 
             dataGenerator = ggmGenerator(ordinal = TRUE, nLevels = 5),
             default = "EBICglasso",
             nCases = c(100,250,500,1000,2500),
             tuning = 0.5,
             nReps = 100,
             nCores = 1
             )

## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).

## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).
## An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).

“# An empty network was selected to be the best fitting network. Possibly set ‘lambda.min.ratio’ higher to search more sparse networks. You can also change the ‘gamma’ parameter to improve sensitivity (at the cost of specificity).”

“## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal = ## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 * ## lambda.max). Recent work indicates a possible drop in specificity. Interpret the ## presence of the smallest edges with care. Setting threshold = TRUE will enforce ## higher specificity, at the cost of sensitivity.”

Sensitivity: Also termed the true-positive rate, the proportion of edges present in the true network that were detected in the estimated network.

Specificity: Also termed the true-negative rate, the proportion of missing edges in the true network that were also detected correctly to be absent edges in the estimated network.

The correlation between edge weights of the true network and edge weights of the estimated network, or between centrality estimates based on the true network and centrality estimates based on the estimated network

trueNetwork$graph

##               1            2           3           4            5          6
## 1   0.000000000  0.195283249 0.005102409  0.27464323  0.093370806 0.00000000
## 2   0.195283249  0.000000000 0.457258512  0.00000000  0.017356395 0.00000000
## 3   0.005102409  0.457258512 0.000000000  0.21631094  0.141626695 0.00000000
## 4   0.274643230  0.000000000 0.216310940  0.00000000  0.303173116 0.03679962
## 5   0.093370806  0.017356395 0.141626695  0.30317312  0.000000000 0.18185716
## 6   0.000000000  0.000000000 0.000000000  0.03679962  0.181857157 0.00000000
## 7   0.000000000  0.067998717 0.015047013  0.12339396  0.122698337 0.22094378
## 8   0.000000000  0.000000000 0.000000000  0.00000000  0.000000000 0.09662244
## 9   0.000000000  0.000000000 0.000000000  0.00000000  0.000000000 0.00000000
## 10  0.000000000  0.000000000 0.000000000  0.09450993 -0.010580349 0.01414188
## 11  0.205659921 -0.010705344 0.021628378  0.07974091  0.007534505 0.00000000
## 12  0.064983356 -0.008753877 0.098501073  0.00000000  0.000000000 0.00000000
## 13  0.000000000  0.000000000 0.000000000  0.00000000  0.000000000 0.00000000
## 14  0.000000000  0.000000000 0.000000000 -0.20729516  0.000000000 0.00000000
## 15 -0.005972199  0.000000000 0.000000000  0.00000000  0.118998463 0.00000000
## 16  0.000000000  0.130339925 0.000000000  0.00000000  0.000000000 0.00000000
## 17  0.000000000  0.000000000 0.076840447  0.00000000  0.000000000 0.05562034
## 18  0.000000000  0.064252350 0.072374885  0.00000000  0.101640680 0.00000000
## 19  0.000000000  0.000000000 0.000000000  0.00000000  0.000000000 0.00000000
## 20  0.007479995  0.113931755 0.096802652  0.00000000  0.036211153 0.00000000
##             7          8          9          10           11           12
## 1  0.00000000 0.00000000 0.00000000  0.00000000  0.205659921  0.064983356
## 2  0.06799872 0.00000000 0.00000000  0.00000000 -0.010705344 -0.008753877
## 3  0.01504701 0.00000000 0.00000000  0.00000000  0.021628378  0.098501073
## 4  0.12339396 0.00000000 0.00000000  0.09450993  0.079740908  0.000000000
## 5  0.12269834 0.00000000 0.00000000 -0.01058035  0.007534505  0.000000000
## 6  0.22094378 0.09662244 0.00000000  0.01414188  0.000000000  0.000000000
## 7  0.00000000 0.03033601 0.00000000  0.00000000  0.000000000  0.106636372
## 8  0.03033601 0.00000000 0.00000000  0.09123952  0.022549418  0.000000000
## 9  0.00000000 0.00000000 0.00000000  0.00000000  0.271642439  0.146712454
## 10 0.00000000 0.09123952 0.00000000  0.00000000  0.505686959  0.000000000
## 11 0.00000000 0.02254942 0.27164244  0.50568696  0.000000000  0.000000000
## 12 0.10663637 0.00000000 0.14671245  0.00000000  0.000000000  0.000000000
## 13 0.00000000 0.00000000 0.03047257  0.00000000  0.000000000  0.129588198
## 14 0.00000000 0.00000000 0.15097476  0.00000000  0.000000000  0.118422046
## 15 0.08505491 0.02048741 0.07485573  0.00000000  0.180410481  0.067038344
## 16 0.04628658 0.12054616 0.00000000  0.03115945  0.000000000  0.000000000
## 17 0.00000000 0.06895300 0.02888292  0.00000000  0.000000000  0.000000000
## 18 0.00000000 0.00000000 0.00000000  0.00000000  0.000000000  0.010957557
## 19 0.00000000 0.10858793 0.04819525  0.00000000  0.024906667  0.299685751
## 20 0.00000000 0.00000000 0.00000000  0.00000000  0.000000000  0.000000000
##            13          14           15         16           17         18
## 1  0.00000000  0.00000000 -0.005972199 0.00000000  0.000000000 0.00000000
## 2  0.00000000  0.00000000  0.000000000 0.13033993  0.000000000 0.06425235
## 3  0.00000000  0.00000000  0.000000000 0.00000000  0.076840447 0.07237489
## 4  0.00000000 -0.20729516  0.000000000 0.00000000  0.000000000 0.00000000
## 5  0.00000000  0.00000000  0.118998463 0.00000000  0.000000000 0.10164068
## 6  0.00000000  0.00000000  0.000000000 0.00000000  0.055620338 0.00000000
## 7  0.00000000  0.00000000  0.085054913 0.04628658  0.000000000 0.00000000
## 8  0.00000000  0.00000000  0.020487413 0.12054616  0.068953005 0.00000000
## 9  0.03047257  0.15097476  0.074855727 0.00000000  0.028882922 0.00000000
## 10 0.00000000  0.00000000  0.000000000 0.03115945  0.000000000 0.00000000
## 11 0.00000000  0.00000000  0.180410481 0.00000000  0.000000000 0.00000000
## 12 0.12958820  0.11842205  0.067038344 0.00000000  0.000000000 0.01095756
## 13 0.00000000  0.50456045  0.011820020 0.06390397  0.089867535 0.09355390
## 14 0.50456045  0.00000000  0.153702903 0.00000000  0.000000000 0.00000000
## 15 0.01182002  0.15370290  0.000000000 0.27661091 -0.002980352 0.00000000
## 16 0.06390397  0.00000000  0.276610912 0.00000000  0.036502063 0.11124936
## 17 0.08986753  0.00000000 -0.002980352 0.03650206  0.000000000 0.34457546
## 18 0.09355390  0.00000000  0.000000000 0.11124936  0.344575459 0.00000000
## 19 0.15698509  0.02847355 -0.014459351 0.00000000  0.000000000 0.08224302
## 20 0.10422255  0.04081857  0.000000000 0.00000000  0.007031625 0.10910560
##             19          20
## 1   0.00000000 0.007479995
## 2   0.00000000 0.113931755
## 3   0.00000000 0.096802652
## 4   0.00000000 0.000000000
## 5   0.00000000 0.036211153
## 6   0.00000000 0.000000000
## 7   0.00000000 0.000000000
## 8   0.10858793 0.000000000
## 9   0.04819525 0.000000000
## 10  0.00000000 0.000000000
## 11  0.02490667 0.000000000
## 12  0.29968575 0.000000000
## 13  0.15698509 0.104222550
## 14  0.02847355 0.040818573
## 15 -0.01445935 0.000000000
## 16  0.00000000 0.000000000
## 17  0.00000000 0.007031625
## 18  0.08224302 0.109105596
## 19  0.00000000 0.201752964
## 20  0.20175296 0.000000000

head(simRes)

## Warning: `select_()` is deprecated as of dplyr 0.7.0.
## Please use `select()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## Warning: `summarise_each_()` is deprecated as of dplyr 0.7.0.
## Please use `across()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## === netSimulator Results ===
## 
## Mean (SD) values per varied levels:
## 
##   nCases    default tuning ExpectedInfluence MaxFalseEdgeWidth       bias
## 1    100 EBICglasso    0.5         0.8500182         0.3011283 0.04014474
## 2    500 EBICglasso    0.5         0.9213254         0.2011513 0.02145479
## 3   1000 EBICglasso    0.5         0.9549089         0.1258747 0.01717972
## 4   1000 EBICglasso    0.5         0.9714165         0.1874430 0.01843795
## 5   1000 EBICglasso    0.5         0.9763828         0.1581449 0.01749627
## 6   2500 EBICglasso    0.5         0.9818661         0.1106360 0.01429199
##   sensitivity specificity correlation  strength closeness betweenness
## 1   0.62 (NA)   0.72 (NA)   0.71 (NA) 0.81 (NA) 0.42 (NA)   0.49 (NA)
## 2   0.89 (NA)   0.62 (NA)   0.94 (NA) 0.93 (NA) 0.85 (NA)   0.92 (NA)
## 3   0.86 (NA)   0.68 (NA)   0.96 (NA) 0.95 (NA) 0.84 (NA)   0.87 (NA)
## 4   0.89 (NA)   0.64 (NA)   0.96 (NA) 0.98 (NA) 0.87 (NA)   0.86 (NA)
## 5   0.92 (NA)   0.65 (NA)   0.96 (NA) 0.97 (NA) 0.94 (NA)   0.93 (NA)
## 6   0.94 (NA)   0.51 (NA)   0.98 (NA) 0.99 (NA) 0.94 (NA)    0.8 (NA)
## 
## 
## Use plot(x) to plot results (nCases only), or as.data.frame(x) to see all results.

## Figure:
#pdf("Fig6.pdf",width=8,height=4)
plot(simRes)

# N=250 achieves a correlation between the “true” and estimated networks above 0.8 for edge weights and strength, and above 0.7 for sensitivity.


plot(simRes,yvar = c("strength","closeness","betweenness"))

#dev.off()

"We recommend:

1. more research estimating network models from psychological data will make clear what one could expect as a true network structure, especially if researchers make the statistical parameters of their network models publicly available;

2. researchers should simulate network models under a wide variety of potential true network structures, using different estimation methods;

3. researchers should simulate data under an expected network structure to gain some insight in the required sample size. "

Post-hoc Stability checks

ncores <- 1
if(detectCores()-1 > 1) {
    ncores <- detectCores()-5 # modify
}


boot1 <- try(bootnet(results, nCores = ncores, nBoots = 1000, type = "nonparametric"))

## Note: bootnet will store only the following statistics:  edge, strength, outStrength, inStrength

## Error in is(data, "bootnetResult") : 找不到对象'results'

boot1 <- try(bootnet(data, nCores = ncores, nBoots = 1000, type = "nonparametric"))

## Note: bootnet will store only the following statistics:  edge, strength, outStrength, inStrength

## Estimating sample network...

## Warning in formals(fun): argument is not a function

## Warning in formals(fun): argument is not a function

## Error in do.call(.input$estimator, c(list(data), .input$arguments)) : 
##   'what' must be a function or character string

graph3 <- estimateNetwork(
  data,
  default = "EBICglasso",
  corMethod = "cor_auto",
  tuning = 0.5)

boot1 <- bootnet(graph3, nCores = ncores, nBoots = 1000, type = "nonparametric")
boot2 <- bootnet(graph3, nCores = ncores, nBoots = 1000, type = "case", statistics=c("strength", "closeness","betweenness"))

# nonparametric bootstrap (using resampled data with replacement)
# Confidence intervals can not be constructed for centrality indices (see the supplementary materials of Epskamp, Borsboom et al., 2017). 
# To assess the stability of centrality indices, one can perform a case-dropping bootstrap (subsampling without replacement).

## Bootstrap plots:
#pdf("BootnetResults.pdf", height = 10, width = 15)
plot(boot1, order = "sample")

## Warning: `filter_()` is deprecated as of dplyr 0.7.0.
## Please use `filter()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## Warning: `summarise_()` is deprecated as of dplyr 0.7.0.
## Please use `summarise()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## Warning: `mutate_()` is deprecated as of dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

plot(boot1, "edge", plot = "difference", onlyNonZero = TRUE,
     order = "sample")

## Expected significance level given number of bootstrap samples is approximately: 0.05

plot(boot1, "strength", order = "sample", plot = "difference")

## Expected significance level given number of bootstrap samples is approximately: 0.051

plot(boot2)

## CS-coefficient:
corStability(boot2)

## === Correlation Stability Analysis === 
## 
## Sampling levels tested:
##    nPerson Drop%   n
## 1       55  75.1  17
## 2       72  67.4  82
## 3       90  59.3  95
## 4      107  51.6 135
## 5      124  43.9 130
## 6      141  36.2 116
## 7      158  28.5 103
## 8      176  20.4 128
## 9      193  12.7  96
## 10     210   5.0  98
## 
## Maximum drop proportions to retain correlation of 0.7 in at least 95% of the samples:
## 
## betweenness: 0.05 (CS-coefficient is lowest level tested)
##   - For more accuracy, run bootnet(..., caseMin = 0, caseMax = 0.127) 
## 
## closeness: 0.127 
##   - For more accuracy, run bootnet(..., caseMin = 0.05, caseMax = 0.204) 
## 
## strength: 0.516 
##   - For more accuracy, run bootnet(..., caseMin = 0.439, caseMax = 0.593) 
## 
## Accuracy can also be increased by increasing both 'nBoots' and 'caseN'.